First, we load the data into the workspace.
Following is a summary of the dataset.
summary(checkin.NY)
## weather.id gid user_id gender
## Min. : 24 Min. : 1 Min. : 32 female:259667
## 1st Qu.: 888 1st Qu.: 276193 1st Qu.: 1542768 male :311595
## Median :1763 Median : 552188 Median : 6913375 NA's : 8524
## Mean :1760 Mean : 551507 Mean :15589057
## 3rd Qu.:2589 3rd Qu.: 825534 3rd Qu.:21767105
## Max. :3590 Max. :1135513 Max. :89322077
##
## venue_id lat lon
## 43a52546f964a520532c1fe3: 12663 Min. :40.5 Min. :-74.2
## 4ace6c89f964a52078d020e3: 7886 1st Qu.:40.7 1st Qu.:-74.0
## 42911d00f964a520f5231fe3: 3521 Median :40.7 Median :-74.0
## 4ae6363ef964a520aba521e3: 2967 Mean :40.7 Mean :-74.0
## 49b7ed6df964a52030531fe3: 2835 3rd Qu.:40.8 3rd Qu.:-74.0
## 42829c80f964a5206a221fe3: 2191 Max. :40.9 Max. :-73.7
## (Other) :547723
## timestamps.x cate_l1
## Min. :1.39e+09 Food :170208
## 1st Qu.:1.39e+09 Shop & Service : 90591
## Median :1.40e+09 Nightlife Spot : 70389
## Mean :1.40e+09 Travel & Transport : 62418
## 3rd Qu.:1.40e+09 Arts & Entertainment : 60727
## Max. :1.40e+09 Professional & Other Places: 57076
## (Other) : 68377
## cate_l2 datetime
## Gym / Fitness Center: 32386 Min. :2014-02-01 00:00:04
## Bar : 25039 1st Qu.:2014-03-09 19:27:56
## Airport : 24891 Median :2014-04-15 11:12:26
## Office : 24356 Mean :2014-04-15 06:30:21
## Coffee Shop : 18346 3rd Qu.:2014-05-20 03:09:42
## American Restaurant : 15471 Max. :2014-06-30 23:58:52
## (Other) :439297 NA's :69
## hour yearday weekday isweekend
## 19 : 57300 Length:579786 Length:579786 Weekend:192831
## 18 : 52049 Class :character Class :character Workday:386955
## 20 : 45170 Mode :character Mode :character
## 13 : 38504
## 17 : 37898
## 14 : 35523
## (Other):313342
## last.gid last.timestamps last.venue_id
## Min. : -1 Min. :-1.00e+00 43a52546f964a520532c1fe3: 9083
## 1st Qu.: 203258 1st Qu.: 1.39e+09 4ace6c89f964a52078d020e3: 6081
## Median : 491981 Median : 1.40e+09 42911d00f964a520f5231fe3: 2972
## Mean : 498473 Mean : 1.29e+09 4ae6363ef964a520aba521e3: 2628
## 3rd Qu.: 784803 3rd Qu.: 1.40e+09 49b7ed6df964a52030531fe3: 2456
## Max. :1135513 Max. : 1.40e+09 (Other) :513535
## NA's : 43031
## last.lat last.lon last.cate_l1
## Min. :-1.0 Min. :-74.2 Food :158896
## 1st Qu.:40.7 1st Qu.:-74.0 Shop & Service : 86910
## Median :40.7 Median :-74.0 Nightlife Spot : 65845
## Mean :37.6 Mean :-68.5 Arts & Entertainment : 54447
## 3rd Qu.:40.8 3rd Qu.:-73.9 Professional & Other Places: 53443
## Max. :40.9 Max. : -1.0 (Other) :117214
## NA's : 43031
## last.cate_l2 last.user_id time.interval
## Gym / Fitness Center: 31436 Min. : 32 Min. : 0
## Bar : 23525 1st Qu.: 1541303 1st Qu.: 3
## Office : 22786 Median : 6812580 Median : 21
## Airport : 18457 Mean :15394230 Mean : 28852
## Coffee Shop : 17437 3rd Qu.:21345141 3rd Qu.: 93
## (Other) :423114 Max. :89322077 Max. :390052
## NA's : 43031 NA's :43031
## latitude longitude conds winddird
## Min. :40.7 Min. :-73.9 Clear :289609 Min. : 0.0
## 1st Qu.:40.7 1st Qu.:-73.9 Overcast :134203 1st Qu.: 0.0
## Median :40.7 Median :-73.9 Mostly Cloudy : 46510 Median : 0.0
## Mean :40.7 Mean :-73.9 Partly Cloudy : 31744 Mean : 92.7
## 3rd Qu.:40.7 3rd Qu.:-73.9 Scattered Clouds: 28311 3rd Qu.:190.0
## Max. :40.7 Max. :-73.9 (Other) : 47843 Max. :360.0
## NA's : 1566
## windspd temperatur fog rain
## Min. : 0 Min. :-12.2 Mode :logical Mode :logical
## 1st Qu.: 6 1st Qu.: 4.4 FALSE:576316 FALSE:551962
## Median : 9 Median : 12.8 TRUE :3470 TRUE :27824
## Mean :10 Mean : 11.8 NA's :0 NA's :0
## 3rd Qu.:13 3rd Qu.: 19.4
## Max. :35 Max. : 31.1
## NA's :38929 NA's :188
## snow thunder tornado
## Mode :logical Mode :logical Mode :logical
## FALSE:573449 FALSE:579786 FALSE:579786
## TRUE :6337 NA's :0 NA's :0
## NA's :0
##
##
##
freq.plot(checkin.NY,title="New York City")
The analysis in the frequency domain reveals a strong harmonics with a period of 24 hours. Therefore, it is reasonable to analyze the data by 24 hours.
The distribution of the check-in categories accross one period (24 hours).
time.distribution.plot(checkin.NY,title="New York City")
The distribution of the check-in categories in 24 hours in radial plots.
time.radial.plot(checkin.NY)
have a look at a intuitive spatiotemporal visualizaiton of the data…
basemap = map.plot("../../global/data/shapefiles",
"NYC_borough_boundaries_WGS84",
alpha=0.1,size=0.3,color="grey")
saveGIF({
point.animation.plot(checkin.NY,title="New York City",
color="cate_l1", color.concrete=FALSE,
basemap=basemap,axis.size=18,
title.size=24,legend.size=18,point.size=3)
}, interval = 0.5, movie.name = "spatiotemporal.gif",
ani.width = 1500, ani.height = 1200)
What will the statistics say when we neglect the temporal factors? First of all, if we try to find spatial clusters based on differnt checkin categories, the distribution of the founded clusters are quite differnt. It indicates that each category has differnt correlations with the geographic space.
if(!exists("basemap")){
basemap <- map.plot("../../global/data/shapefiles",
"NYC_borough_boundaries_WGS84",
alpha=0.1,size=0.3,color="grey")
}
## OGR data source with driver: ESRI Shapefile
## Source: "../../global/data/shapefiles", layer: "NYC_borough_boundaries_WGS84"
## with 5 features and 4 fields
## Feature type: wkbMultiPolygon with 2 dimensions
plots <- lapply(split(checkin.NY,checkin.NY$cate_l1),function(ci){
# find out clusters for the type of category
clusters = spatial.clustering(ci)
centers = clusters[["centers"]]
points = clusters[["point.unique"]]
wss = clusters[["wss"]]
pct = clusters[["pct"]]
lbl = paste(nrow(centers),"clusters;\nWSS:",
formatC(wss,digits=2,format="f"),
"(",format.percent(pct),")")
# add basic points
gg.map <- point.plot(points,x="lon.x",y="lat.x",alpha=0.05, basemap=basemap)
# add cluster information
gg.map <- point.plot(centers, x="lon.center",y="lat.center",
color= "cid.ordered", color.concrete=FALSE,
point.size = log(centers$size,5),alpha = 0.7,
basemap=gg.map,
xlim=range(checkin.NY$lon,finite=TRUE),
ylim=range(checkin.NY$lat,finite=TRUE))
# some plot configuration
gg.map <- gg.map + ggtitle(ci[1,"cate_l1"]) +
theme(legend.position="none") +
geom_text(aes(x = -74.13, y = 40.87), label = lbl, size=2)
})
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
# png("categorized_clustering.png",width=15*ppi,height=6*ppi,res=ppi)
# grid.arrange(plots[[1]],plots[[2]],plots[[3]],plots[[4]],
# plots[[5]],plots[[6]],plots[[7]],plots[[8]],
# plots[[9]],plots[[10]],nrow=4, ncol=3)
# dev.off()
plots[[1]];plots[[2]];plots[[3]];plots[[4]];plots[[5]];plots[[6]];plots[[7]];plots[[8]]; plots[[9]];plots[[10]]
We could also see the most domimant categories in each part of the entire city.
# prepare polygon data
SPDF = readOGR(dsn = "../../global/data/shapefiles", layer = "NYC_zipcode")
## OGR data source with driver: ESRI Shapefile
## Source: "../../global/data/shapefiles", layer: "NYC_zipcode"
## with 236 features and 6 fields
## Feature type: wkbPolygon with 2 dimensions
# intersect the checkin data with the each postal code region
checkin.NY = point.in.poly(checkin.NY, SPDF, copy.attr="POSTAL")
## [1] "Spatial data has been prepared."
## [1] "The information from the polygon has been assigned to the points."
# get the distribution of categories in each postal code region
cate.distr=cate.distr.in.poly(checkin.NY)
## [1] "The distribution of category has been assigned to each polygon."
# and use that information to implement the SPDF
SPDF = poly.with.category(SPDF, cate.distr=cate.distr, poly.attr="POSTAL")
# plot
mapdf=df.from.spdf(SPDF)
# mapdf$density=apply(mapdf,1,function(i){i[i["cate.dom"]]})
# mapdf$density=as.numeric(formatC(mapdf$density,digits=1,format = "f"))
# png("main_category.png",height=6*ppi, width=16*ppi,res=ppi)
# grid.arrange(
map.plot(mapdf=mapdf,
more.aes=aes_string(fill="Category.1st"),
color="grey",size=0.3,alpha=0.9)+
theme_bw()+
xlab("")+ylab("")+
theme(legend.title = element_text(size=8),
legend.text = element_text(size = 8))
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
map.plot(mapdf=mapdf,
more.aes=aes_string(fill="Category.2nd"),
color="grey",size=0.3,alpha=0.7)+
theme_bw()+
xlab("")+ylab("")+
theme(legend.title = element_text(size=8),
legend.text = element_text(size = 8))
## Warning: Bedingung hat Länge > 1 und nur das erste Element wird benutzt
# nrow=1, ncol=2, widths=c(1,1) )
# dev.off()
cate.conds = xtabs(~conds+cate_l1, data=checkin.NY)
#prop.table(cate.conds, 1) # row percentages
#prop.table(cate.conds, 2) # column percentages
fit <- ca(cate.conds)
#print(fit) # basic results
summary(fit) # extended results
##
## Principal inertias (eigenvalues):
##
## dim value % cum% scree plot
## 1 0.003244 52.2 52.2 *************************
## 2 0.002082 33.5 85.8 ****************
## 3 0.000511 8.2 94.0 ****
## 4 0.000147 2.4 96.4 *
## 5 7.6e-050 1.2 97.6 *
## 6 6.7e-050 1.1 98.7
## 7 5.1e-050 0.8 99.5
## 8 2.4e-050 0.4 99.9
## 9 8e-06000 0.1 100.0
## -------- -----
## Total: 0.006210 100.0
##
##
## Rows:
## name mass qlt inr k=1 cor ctr k=2 cor ctr
## 1 | Cler | 501 779 88 | 22 452 76 | -19 327 86 |
## 2 | Fog | 2 634 31 | -265 628 37 | -26 6 1 |
## 3 | Haze | 22 975 178 | -200 798 272 | 94 177 94 |
## 4 | HvyR | 5 917 75 | 208 421 61 | 226 496 111 |
## 5 | HvyS | 0 116 5 | -170 20 0 | 377 96 1 |
## 6 | LgFR | 0 714 29 | -590 712 40 | 28 2 0 |
## 7 | LghR | 33 947 88 | 77 362 61 | 98 586 154 |
## 8 | LghS | 7 246 21 | -62 203 8 | 29 44 3 |
## 9 | MstC | 80 738 38 | -31 319 23 | -35 420 47 |
## 10 | Ovrc | 232 957 107 | -25 217 45 | 46 740 237 |
## 11 | PrtC | 55 592 41 | 38 315 25 | -36 277 34 |
## 12 | Rain | 10 943 67 | 163 632 81 | 114 312 62 |
## 13 | SctC | 49 874 97 | -59 287 53 | -85 587 169 |
## 14 | Snow | 4 850 134 | -421 849 218 | 4 0 0 |
##
## Columns:
## name mass qlt inr k=1 cor ctr k=2 cor ctr
## 1 | ArtE | 105 763 76 | 49 526 76 | 33 237 54 |
## 2 | CllU | 17 677 27 | -75 580 30 | 31 97 8 |
## 3 | Evnt | 5 557 20 | -24 23 1 | -113 534 32 |
## 4 | Food | 294 688 48 | 18 301 28 | -20 388 56 |
## 5 | NghS | 122 955 206 | 93 826 325 | 37 129 80 |
## 6 | OtdR | 69 890 200 | -9 4 2 | -127 886 528 |
## 7 | PrOP | 98 986 272 | -123 878 457 | 43 108 87 |
## 8 | Rsdn | 27 352 42 | 23 54 4 | 54 298 38 |
## 9 | ShpS | 156 661 49 | -29 430 41 | -21 230 34 |
## 10 | TrvT | 107 791 59 | -33 318 36 | 40 473 84 |
#plot(fit) # symmetric map
plot(fit, mass = TRUE, contrib = "absolute", map ="rowgreen",
arrows = c(TRUE, FALSE)) # asymmetric map
Under the assumption \(H=i\) is independent from \(W=j\) \[P(C=k|H=i,W=j)=\frac{P(C=k,H=i,W=j)}{P(H=i,W=j)}=\frac{P(H=i,W=j|C=k)*P(C=k)}{P(H=i)*P(W=j)} (1) \]
since \(H=i\) is independent from \(W=j\), \[Exp[P(H=i,W=j|C=k)]=P(H=i|C=k)*P(W=j|C=k) (2)\]
therefore, \[Exp[P(C=k|H=i,W=j)]=Exp[ \frac{P(H=i,W=j|C=k)*P(C=k)} {P(H=i)*P(W=j)}] \\\ =\frac{P(H=i|C=k)*P(W=j|C=k)*P(C=k)}{P(H=i)*P(W=j)} \\\ =\frac{\frac{P(H=i,C=k)}{P(C=k)}*\frac{P(W=j,C=k)}{P(C=k)}*P(C=k)}{P(H=i)*P(W=j)} \\\ =\frac{P(C=k|H=i)*P(H=i)*P(C=k|W=j)*P(W=j)}{P(H=i)*P(W=j)*P(C=k)} \\\ =\frac{P(C=k|H=i)*P(C=k|W=j)}{P(C=k)} (3)\]
\[P_{u}(C=k|H=i)=\frac{\Phi_{u}(C=k,H=i)}{\Phi_{u} (H=i)} (4)\]
get.temporal.impact <- function(dataframe,hour){
# dataframe.in.hour = checkin.single[which(checkin.single$hour==hour),]
dataframe.in.hour = dataframe[which(dataframe$hour==hour),]
phi.h = nrow(dataframe.in.hour)
list.category = split(dataframe.in.hour, dataframe.in.hour$cate_l2)
sapply(list.category, function(i){
nrow(i)/phi.h
})
}
\[P_{u}(C=k|W=j)=\frac{Intercept(C=k,W=j)}{\sum Intercept(C,W=j)} (5)\]
get.meteorological.impact <- function(fit,conds){
conds.id = which(fit[["rownames"]]==conds)
ref.vec = fit[["rowcoord"]][conds.id,1:8]
cate.all = fit[["colcoord"]][,1:8]
intercepts = apply(cate.all, 1, function(x){
(x[1]*ref.vec[1] + x[2]*ref.vec[2] + x[3]*ref.vec[3] + x[4]*ref.vec[4] +
x[5]*ref.vec[5] + x[6]*ref.vec[6] + x[7]*ref.vec[7] + x[8]*ref.vec[8] ) /
(ref.vec[1]^2 + ref.vec[2]^2 + ref.vec[3]^2 + ref.vec[4]^2 +
ref.vec[5]^2 + ref.vec[6]^2 + ref.vec[7]^2 + ref.vec[8]^2 )
} )
# vec = intercepts / sum(intercepts) !!!wrong!!! intercetp can be negative
vec = (intercepts - min(intercepts)) / (max(intercepts)-min(intercepts)) # scale into [0,1]
names(vec) = fit[["colnames"]]
vec
}
get.meteorological.impact2 <- function(dataframe,conds){
dataframe.in.conds = dataframe[which(dataframe$conds==conds),]
phi.c = nrow(dataframe.in.conds)
list.category = split(dataframe.in.conds, dataframe.in.conds$cate_l2)
sapply(list.category, function(i){
nrow(i)/phi.c
})
}
\[P_{u}^{*} (C=k|W=j)= w_{j}*[P_{u}(C=k|W=j)-\bar P_{u}]+\bar P_{u}\]
get.weather.weight <- function(fit){
conds.all = fit[["rowcoord"]][,1:2]
mag = apply(conds.all, 1, function(x){
sqrt( (x[1]^2+x[2]^2) )
} )
mag / sum(mag)
}
get.weighted.meteorological.impact <- function(fit, conds, weight){
# cate.conds = xtabs(~conds+cate_l2, data=dataframe)
# fit <- ca(cate.conds)
unweighted = get.meteorological.impact(fit,conds)
# weights = get.weather.weight(fit)
# conds.id = which(fit[["rownames"]]==conds)
# vec = weights[conds.id] * (unweighted-mean(unweighted)) + mean(unweighted)
vec = weight * (unweighted-mean(unweighted)) + mean(unweighted)
names(vec) = fit[["colnames"]]
vec
}
get.weighted.meteorological.impact2 <- function(dataframe, conds, weight){
unweighted = get.meteorological.impact2(dataframe, conds)
weight * (unweighted-mean(unweighted)) + mean(unweighted)
}
\[P(C=k)=\frac{\Phi_{u} (C=k) }{\Phi_{u} }\]
get.denominator <- function(dataframe){
phi.h = nrow(dataframe)
list.category = split(dataframe, dataframe$cate_l2)
denominator = sapply(list.category, function(i){
nrow(i)/phi.h
})
denominator
}
\[E[P(C=k|H=i,W=j)]=\frac{P(C=k|H=i)*P(C=k|W=j)}{P(C=k)}\]
update @2014.09.04
p.k has a relatively high performance, multipling it with others do not bringer a higher one; instead, dividing does.update @2014.09.05
get.meteorological.impact(), I used intercept/sum(intercept), where intercept can be negative, and the sum(intercept) can be nonsense. Now it is fixed to (intercepts - min(intercepts)) / (max(intercepts)-min(intercepts)) to scale into [0,1].p.kj is calculated by just using the probability distribution, then the equation p.kij = p.ki * p.kj / p.k is correct (which confirmes our derivation);p.kj is calculated by personal CA, it should be noticed that the scaled interception does not describe the probability of C=k under the condition W=j. Instead, it describes the negative/positive impacts. Hence, the conditional probability should be the scaled interception multipled by the priori p.k. Therefore, here the final result should be p.kij = p.ki * p.kj.p.j, the above two ways have very similar results (which is expected!).